Momentum Average approximation for models with boson-modulated hopping: the 
role of closed loops in the dynamical generation of a finite quasiparticle mass 



Mona Berciu 1 and Holger Fehske 2 

1 Department of Physics and Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada 
2 Institut fur Physik, Ernst-Moritz-Arndt-Universitat Greifswald, D-17487 Greifswald, Germany 

(Dated: October 21, 2010) 

We generalize the momentum average (MA) approximation to study the properties of single po- 
larons in models with boson affected hopping, where the fermion-boson scattering depends explicitly 
on both the fermion's and the boson's momentum. As a specific example, we investigate the Edwards 
fermion-boson model in both one- and two-dimensions. In one dimension, this allows us to compare 
our results with Exact Diagonalization results, to validate the accuracy of our approximation. The 
generalization to two-dimensional lattices allows us to calculate the polaron's quasiparticle weight 
and dispersion throughout the Brillouin zone, and to demonstrate the importance of Trugman loops 
in generating a finite effective mass even when the free fermion has an infinite mass. 
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I. INTRODUCTION 

One of the most common problems in condensed mat- 
ter physics is that of understanding the behavior of a 
particle coupled to bosons from its environment, for ex- 
ample an electron interacting with phonons, magnons, 
or orbitons of its host crystal^— The particle becomes 
"dressed" by a cloud of bosonic excitations that accom- 
pany it. The resulting composite object, generally known 
as a polaron, can have properties significantly different 
from that of the bare particle^— 

The theoretical study of these properties is rather diffi- 
cult away from the various asymptotic regimes where per- 
turbation theory holds. Of course, a large variety of nu- 
merical techniques have been developed to deal with such 
problems jiSr— and many interesting results have been 
uncovered, although the focus so far has been primarily 
on rather simple models such as the Holstein Hamilto- 
nia n 17 : 18 that describes the simplest possible electron- 
phonon coupling. The progress on analytical approxima- 
tions that can efficiently yet accurately describe the non- 
perturbative regimes has been slower. In fact, it is only 
recently that the so-called Momentum Average (MA) ap- 
proximation has been proposed for the Holstein model, 
and shown to accurately capture its polaronic behavior 
in all the parameter space except the extreme adiabatic 
limit ^ A way to systematically improve this approxima- 
tion, as well as generalizations to certain kinds of more 
complex models have been proposed since<22r— The avail- 
ability of such simple yet accurate approximations is im- 
portant, as it allows one to quickly explore large regions 
of the parameter space to identify the interesting prop- 
erties of the model. 

In this work we present the generalization of MA-type 
methods to calculate single polaron Green's functions 
for Hamiltonians whose hopping is boson affected. The 
bosons are assumed to be dispersionless, i.e. of Einstein 
type. For most polaron models, including the one dis- 
cussed here, the spin of the fermion is irrelevant and 



we ignore it. Exceptions occur, for example, in sys- 
tems with spin-orbit coupling, where suitable general- 
izations can be implemented^ The fermion moves on a 
d-dimensional lattice, which for simplicity is assumed to 
be hypercubic (generalization to other types of lattices 
is straightforward 23 ). The cases d — 1 and d = 2 for 
the Edwards fermion-boson model^l are discussed in de- 
tail and interesting physics related to the role of closed 
loops, possible in two dimensions (2D) but absent in ID, 
is uncovered. We note that single polaron properties for 
this model have been investigated numerically in lDj 2 ^ 
and we use these results to assess the accuracy of MA. 
We then extend our method to 2D, where no results are 
currently available, and where we illustrate interesting ef- 
fects of the boson modulated hopping. Other such models 
can be treated similarly. 

The general form of the Hamiltonian of interest is: 

k q 

k,q 

Here, Ck and 6 q are fermion, respectively boson anni- 
hilation operators, and N is the number of sites in the 
system. In all our results we assume periodic boundary 
conditions and let N oo, however finite size systems 
and/or other types of boundary conditions can be treated 
similarly. The sums are over the Brillouin zone, ek is the 
free-fermion dispersion while f2 is the bosons' energy (we 
set Ti = 1). Note that in ([1]) the fermion-boson scat- 
tering depends explicitly on both the fermion's and the 
boson's momentum. This is to be contrasted with sim- 
pler cases, such as the Holstein model, where <?(k, q) — > g 
is a constant, or models where the bosons modulate only 
on-site energies but not the hopping integrals, for which 
g(k, q) — > g(q). The accuracy of MA approximations 
for these simpler ty pes of Hamiltonians has been demon- 
strated in Refs. [TOpM 
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We are interested in calculating the single polaron 
Green's function, denned as 

G(k,^) = (0|c k GH4|0), (2) 

where |0) is the vacuum, G(u) = [lu — T~L + ii]}^ 1 is the 
resolvent associated with the Hamiltonian T-L 1 and r\ is a 
positive, infinitesimally small number. 

From the Green's function we get the spectral weight 

A(k,w) = --ImG(k,w) (3) 

which is measurable by (Inverse) Angle-Resolved Pho- 
toemission Spectroscopy^ The lowest-energy pole of 
A(k, lu) allows us to identify the polaron dispersion E(k). 
Its residue is the quasiparticle (qp) weight, 

^ k =|(^| C ] c |0)| 2 , (4) 

i.e. the overlap between the polaron eigenfunction |0k), 
where "H|</>k) = -E(k)|0k)> and a free-fermion state. Of 
course, the spectral weight contains information about 
the higher energy states as well, but here we will primar- 
ily focus on the low-energy polaron band. 

The article is organized as follows. We first introduce 
the Edwards fermion-boson model, which is the specific 
model with boson-modulated hopping that we will use as 
an example in this work. We then outline the MA ap- 
proximation for the simpler ID case, and use comparison 
with available numerical results to analyze its accuracy in 
various regimes. Then, we generalize MA for the 2D case 
and use it to understand the relevance of closed loops for 
generating a dynamical mass for the dresses quasiparti- 
cle. Finally we summarize our results and conclude. 

II. THE EDWARDS FERMION-BOSON MODEL 

The Edwards fermion-boson mode l 27 ' 28 is defined by 
the Hamiltonian 

u = -*f E c h + n E b l b i ~ *«■ E c h ($ + h) , (5) 

where the first term describes nearest-neighbor (NN) 
hopping of the fermion on the lattice of interest, the sec- 
ond term describes the Einstein boson branch, and the 
last term is the boson-modulated hopping. Note that in 
the limit tf — > 0, it is only the last term that allows 
the fermion to move: hopping from one site to a neigh- 
boring one either creates an excitation at the "depar- 
ture" site, or removes one from the "arrival" site. This 
model provides a way to mimic, for example, the mo- 
tion of a fermion through an antiferromagnetically or- 
dered spin background ] 3 ' 30 ' 33 For a Neel antiferromagnet 
(AFM) doped with one fermion — which is the zero-order 
description of a hole moving in a cuprate Cu02 layer — 
the hopping of the fermion reshuffles the spins along its 



(a) 




FIG. 1. Sketch of the 3-boson, 3-site sequence of processes 
that give rise to effective second NN fermion hopping. The 
site occupied by the fermion is shaded and the arrow indicates 
the direction of the next ^-hopping process. Each ^-hopping 
either leaves a boson (drawn as a square) at the initial site, 
or absorbs a boson from the arrival site. Both collinear (a) 
and closed loop (b) processes are allowed in 2D; in ID only 
(a) is possible. 



path. If the spin at the "arrival" site has the proper 
orientation, when it is shuffled to the starting site as 
the fermion hops it will have the wrong orientation (a 
"magnon" defect is created at the initial site). Vice versa, 
if the visited spin has the wrong orientation, when shuf- 
fled by one site it will be properly aligned ("magnon" 
defect removed from the arrival site). This is precisely 
the type of physics described by this boson-modulated 
hopping, although it ignores details such as the hard- 
core boson constraint for magnons, the fact that in a 
Neel AFM the energy of neighboring magnons is not addi- 
tive, and also it allows the particle to coexist with bosons 
at the same site. The free-fermion hopping term tf, on 
the other hand, has to be added when describing motion 
through a Heisenberg-like AFM, where spin fluctuations 
continuously create and annihilate magnons. Indeed, as 
shown in Refs. [28| and [34], the free-fermion hopping term 
can be mapped into a purely bosonic term of the type 
A Ylifii + frj) with A = tfil/ (2i{,), which allows the num- 
ber of magnons to fluctuate. 

The "conventional wisdom" is that a fermion in a 2D 
Neel AFM (tf = case) has an infinite effective mass, 
because as it tries to move it creates a costly string of 
defects which effectively pin it to its original site. This 
is, however, not true. The boson-modulated hopping 
gives rise to an effective fermion mass even when the 
bare fermionic mass is infinite (tf = 0). 28 ' 33 For the ID 
chain, this is primarily due to the 3-site, 3-boson pro- 
cesses sketched in Fig. HJa) which result in an effective 
second NN hopping of the fermion (of course, more com- 
plicated processes involving more bosons are also pos- 
sible, but they are energetically more costly). In Ref. 
I30L it was noted that in 2D, these collinear processes — 
which in 2D give rise to effective 3rd NN hopping — are 
supplemented by the closed loop processes sketched in 
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Fig. [ljb), which give rise to effective 2nd NN hopping. 
The importance of such closed loops — known as Trug- 
man paths — for determining the effective quasiparticle 
mass has been emphasized already in Ref. [3fj in the con- 
text of cuprates. In this work, we are the first to explic- 
itly investigate this phenomenology for the 2D Edwards 
fermion-boson model. 

While we focus on this Hamiltonian for the remainder 
of this work, the MA method can be generalized straight- 
forwardly to other boson-modulated hoppings, like that 
appearing in the ID Su-Schrieffer-Heeger model of poly- 
acetylene^ Its phonon-modulated hopping term is pro- 
portional to J2i (4 c l+ i + H - c -) [ b \ + b i~ b l+i ~ h+i] > in- 
here phonons can be both absorbed and created at either 
of the two sites involved in the hopping process. 

III. THE MOMENTUM AVERAGE 
APPROXIMATION 

One can discuss the meaning of the MA approxima- 
tions from several different points of view. One is that 
this is an approximation which sums semi-analytically all 
diagrams contributing to the self-energy, however "expo- 
nentially small" contributions are ignored when calculat- 
ing the expression of each such diagram^ The precise 
meaning of this statement will be clarified below. 

A more useful starting point for our purposes here is 
the variational meaning of MAi 20 i 31 The central idea of 
polaron physics is that the dressed quasiparticle — the 
polaron — consists of the fermion itself and a cloud of 
bosons in its vicinity. MA permits one to systematically 
select what bosonic states to keep in the variational space 
to describe this cloud, and sum their contributions effi- 
ciently. Of course, the more states kept, the more accu- 
rate the results. Physical intuition is needed to decide 
what is a minimal acceptable starting point. 

In this section, we first describe the MA approximation 
for a ID chain, and assess its accuracy against available 
numerical data. We then briefly review the 2D general- 
ization and analyze the resulting physics, which has not 
been investigated before. 

A. MA for the ID chain 

Since we primarily focus on the low-energy polaronic 
physics, we will only attempt to describe the polaronic 
cloud. The only restriction we impose regards its spatial 
size, i.e., what is the maximum number of neighboring 
sites that it can span. The number of bosons at any site 
in the cloud, as well as the position of the fermion with 
respect to the center of the cloud, are not restricted — 
they can take any values. A few of the possible states for 
a 3-site cloud calculation are sketched in Fig. HJa). Of 
course, any 1- and 2-site cloud states belongs to this set. 

Note that here we do not allow any bosonic excitations 
to occur far from the polaron cloud, since they primarily 




FIG. 2. (Color online) (a) Sketch of several states included in 
the MA 1 - - 1 approximation discussed here. The polaron cloud 
is allowed to extend on up to three neighboring sites, located 
anywhere in the system with any number of bosons (shown as 
red squares) at each site, and the fermion (shown as a shaded 
blue circle) at any distance away from the cloud, (b) Example 
of extra states included in a MA 1 - 1 - 1 level approximation, with 
a boson arbitrarily far away from the polaron cloud. 

contribute to higher-energy states. In other words, this 
will be a MA^ ) level approximation^ Generalization to 
MAW and higher levels, needed for example to describe 
the polaron+one-boson continuum, is straightforward al- 
though fairly cumbersome. For example, in MA*- 1 ' we 
also include states such as sketched in Fig. H^b), with 
one boson arbitrarily far away from the polaron cloud. 
Similarly MA^ 2 ) allows two bosons at arbitrary locations 
away from the cloud, etc. As detailed below, the physics 
encoded in these higher-level approximations does not 
lead to any qualitative changes for the properties and 
parameter ranges that we are interested in. 

After deciding on the MA^ ^ (henceforth called simply 
MA) level, the only question left is how big should one 
allow the polaron cloud to be. For the Holstein model, 
a one-site cloud already gives a remarkably accurate de- 
scription in any dimension, most everywhere in the pa- 
rameter space except intermediate coupling in the adia- 
batic limit Q/t -> 0^22. 

For the Edwards model, as already discussed, at least 
3-site boson clouds need be considered in order to de- 
scribe the leading processes that result in the dynamical 
generation of a finite effective mass in the limit tf, X — > 0. 
As a result, we will start directly by building a variational 
MA approximation allowing any number of bosons on any 
3 consecutive sites, which can be located at any distance 
from where the fermion is. To achieve this, we introduce 
3 types of generalized Green's functions. The first are 

F n (k,q,u>) =J2e i(k - q)R H0\c k G("H b r\V) > (6) 

i 

where the ket describes a state of total momentum k (as 
required since the Hamiltonian is invariant to transla- 
tions) of which the fermion has a momentum q and the 
cloud of n bosons, all located at the same site, has mo- 
mentum k — q. 

Next are the two-site cloud Green's functions, namely 

F 7hm (k,q,u) = Y,e l{k ~ q)R H0\c k G(u;)clbl 7 \^m , 

i 

(7) 
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which are defined only for n > 2 and 1 < m < n — 1, so 
that they are distinct from the one-site cloud functions 

I 



q,u>) defined above. Finally, we have the three-site 
cloud functions, 



F TO (fc,,, w )=^ e ^)^(0| Cfc GM C t & t" 1& t r »-^ +1 |0) 



(8) 



defined for n > 2 and 1 < m < n — 1,1 < p < n — 1 , m + 
p < n. These restrictions again avoid overlap with the 
functions introduced above. Note, however, that they 
allow bosons to exist only on the two outer sites of the 
three-site cloud when m + p = n; such states are not 
accounted for by F n ^ m (k, q,oj). In principle one can keep 
adding other states, cither with more extended clouds, or 
with bosons far away from the cloud, until convergence is 
achieved. For reasons already explained, for this model 
we expect that it suffices to stop here. 

The next step is to generate equations of motion 
for these generalized Green's functions. These are ob- 
tained by using the Dyson identity G(uj) = Gq(uj) + 
G(w)VGo(w) , where Go(w) is the resolvent for Ho = 
H — V. This is an exact (non-perturbative) identity and 

I 



I 

holds for any partitioning of H = Ho + V. For our pur- 
poses, it is convenient to take Ho as the non-interacting 
part and V as the boson-modulated hopping term. Note 
that the kets at the right of G(ui) in all the above defini- 
tions are eigenstates of Ho, so that, for e.g. 

6o(«)cJ&r^T|o> = G («,« - nn)ct 6 r^T|o> , 

(9) 

where 



G (k,Lu) = 



LO + ir]-ek 



(10) 



is the free fermion propagator. With this observation, we 
find the exact equation of motion for G(k, uj) to be 



G(k,uj) = G (k,to) 



ikR, 



(0\c k G(u)(cl 



+i 



)b\\0) 



because Vcf |0) = — 4(cj-i + cj +1 )foj|0). This shows that the Green's function of interest to us is linked to various 
averages over the Brillouin zone (momentum averages) of the F\{k, q, to) Green's functions. To make this more precise, 
we introduce the real-space counterparts of the generalized Green's functions defined above, namely: 



/n(M,w) = -Xy9 Sa F„(fc,<7,") = E 



..ikH, 



(0\c k d{w)cl s b^\0), 



(11) 



^ikRi 



f n , m {k,5,uj) = -J2e^ a F n , m (k,q,cj) =Y,-j W ^\ckG{uj)c\_ 5 b\ &t^™| ) 



(12) 



and 



J kit, 



f m (k,8,Lo) = -Y^e^ a F m (k, q ,u:) = ]T ^(0\c k d(w)c\_ s b^-ibT m ~ P b%i\0) 



(13) 



The kets on the right of the resolvent continue to describe 
a state of total momentum k, however with the fermion 
fixed at a certain distance 8a away from the boson cloud. 
Here a is the lattice distance, so that 5 are integers. These 
definitions do not necessarily lead to the most "esthetic" 
final equations, but they are handy and generalize easily 
to higher dimensions. More symmetric ID equations can 
be found if we use the sin and cos, instead of complex 



r 



Fourier transforms. 

With these notations, we have the exact equation 

G(k, u) = Go(k,iu) [1 - t b (/i(fc, 1, cu) + h{k, -1, u))] . 

(14) 

Let us consider now the equation of motion for 
F n (k,q,uj) with n > 1. When acting on cj6- |0), both 
the boson annihilation and the boson creation part of V 
will contribute. The boson annihilation contribution can 
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be calculated exactly, since bosons can only be annihi- 
lated at the one site where they are present. However, 
the boson creation part can add a boson either at the site 
where the cloud is, or to any other site because a fermion 
in the state cj|0) is delocalized over the entire chain. Of 
course, we do not expect all these outcomes to be equally 
likely, and it is this that allows us to make progress. 



This is where the MA approximation is made: we only 
allow new bosons to be created within at most two-site 
distance from where the one-site boson cloud is. In other 
words, we allow the boson cloud to extend spatially but 
not over more than 3 consecutive sites, since we decided 
that this is our variational space. After all possible such 
terms are accounted for, we find that within this MA 
approximation: 



F n {k, q, u) = -t b G {q, u> - nO) [{e lqa + e'^nf^k, 0, w) + [f n +i(k, 1, w) + f n +i(k, -1, w)] 
+ e~ iqa+ika [/ n+M (k, l,u) + fn+i.i (k, -l,w)] + e lqa [f n +i.n (k,0,u) + /„+i.„(fc, -2,u)] 

+ e- 2l « a+ika {f n+1 sJk,2,u J ) + f n+hhn (k,0,uj)}+e . (15) 



The terms on the first line come from the exact boson 
annihilation part (the first) and the boson creation term 
where the new boson is added at the site where the cloud 
is (the last). The terms in the second line describe the 
contributions where the new boson is created on a NN 



site of the cloud. The terms on the 3rd line are the 
contributions when the new boson is created on a second 
NN site of the cloud. 

Within the same approximation, we also find 



-fn,m {k 



'„q,u) = -t b G (q,u- nil) [(e lqa + e- iqa )mf n -i, m -i(k,0,u) + (1 + e 2lqa )(n - m)f n -i, m {k, -1, u) 
- e iqa [f n+ i, m (k, 0,cd) + / n+ i, m (A;,-2,£j)] + [f n +i, m +i(k, 1,cj) + f n +i, m +i(k, -1, w)] 

" e 2iqa - ika [f n+ l.rn,l(k, 0, u) + /n+l, m ,l(*. " 2 > ^ + ^ [U+l, l,n-m [k, 0, w) + /„ + l,l,„- m (fc, 2, W )]] 



(16) 



and 



F n , m , p (k,q,u;) = -t b G {q,u- nil) [{e~ 2lqa + l)m/„_ 1|m _i, p (fe, l,w) + (e l?a + e' l9Q )(n - m-p)f n -i, miP (k,0,u) 
+ (1 + e 2iqa ) v j n _ lim . v ^ x (k, -1, w) + e- i3 °[/„ +1>ro+1)P (A;, 0, w) + /„+i, m+ i, p (fc, 2, w)] 

+ [/n+l,m,p(fe, 1, W) + fn+l,m,p(k, -1,L))] + e iqa [f n+ i.m iP+ i(k, 0, w) + /n+l, m ,p+l(fc, -2,w)]] . (17) 

As before, creation processes are only allowed to add extra bosons so that the total cloud does not extend over more 
than 3 consecutive sites. The annihilation processes are treated exactly, however one has to be careful when there is 
a single boson on an outside site of the cloud. If this is the annihilated boson, the size of the cloud decreases. As a 
result: 

I 

with a total of n bosons to Green's functions with n — 1 
and n + 1 bosons. The only approximation is the restric- 
tion on the size of the allowed boson cloud. A solution 
of this system will give G(k,u) within this MA approx- 
imation, together with all the other generalized Green's 
functions from which we can extract additional informa- 
tion on the structure of the polaron cloud^ 

This solution is straightforward to obtain in terms of 
the momentum averaged Green's functions f n . First, 
note that only a finite number of these are needed, for a 
given value of n, in order to be able to calculate every- 
thing else. For example, only f n ,m,p(k,S,Lj) with \5\ < 2 
appear on the right hand side of all these equations, and 
similar bounds can be found for the other functions. We 



/„_i,o(M,w) -> e- ife %-i(M + , (18) 

fn-l,n-l(k,S,Lj) -> f n -l(k,S,Cj) , (19) 

fn-i,o, P (k,5,u) -> / n _i.„_i_ p (fc, S,oj), if p< n - 1, 

(20) 

/n-l,0,n-l(M>w) -> e- <fco /n-l(M + l,w) , (21) 

f n -i,mfi{k,5,uj) -> e lka f n ^i. m (k,8 - l,w), if m < n - 1, 

(22) 

/n-i,n-i,o(M,w) -> e ika f n -i(k,6- l,w) . (23) 
All these identities follow directly from the definitions of 

Eqs. nnn-ci- 

We have thus generated an infinite system of coupled 
equations of motion linking various Green's functions 
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therefore first generate a set of recurrence equations for these quantities, using Eqs. (TTTT) - (|T3|) . These read: 
I 

fn{8) = -h [go{8 + l,w„) +g {5 - l,w„)] n/ n _i(0) - i&sro(<*, w„)[/ n+ i(l) + / n +i(-l)] 

- t b g (5 - l,u> n )e ika [/„+!,! (1) + /„+i,i(-l)] - i fc .go(<5 + 1, w„)[/„+i,„(0) + /„+i,„(-2)] 

- t b g (S - 2,Lj n )e +ika [f n+1 ,i, n {2) + /„+i.i,„(0)] - t b g (S + 2,oj n )e- %ka [f n+hn , 1 (0) + /„+i,„,i(-2)] , (24) 



fn, m (S) = -t b [.9o(<5 + l,w„) + ffo(^ - l,w„)]m/„_i, m _i(0) - i b [5o(^w„) +5o(^ + 2,w n )](n - m)/„_i, m (-l) 
— ibffo(<5+ l)W n )[/ n +l m (Q) + /n+l,m( — 2)] — tbgo{5, w n )[/ n +i im +i(l) + / n +i,m+i( — 1)] 



4go(<5 + 2,w„)e- l ' ca [/ 5 



(0) + /n+l,l,n(-2)] ~ *6So(<J - l,Wn)[/i 



n+l,l,n 



-m(0) + /, 



n+l,l,n 



-m(2)] , (25) 



and 



fn, m , P ( s ) = -h\go(S - 2,w„) + 3o(^w n )]m/ n _i >m _i iP (l) - tb[.9o(<5 - l,w„) + .g (<5 + l,o;„)](n - m - p)f n -i, m , P (0) 

- tb{go(5,uj n ) + go(S + 2,w„)]p/ n _i, m j,_i(-l) - t b go(5 - l,u n )[fn+i,m+i,p(0) + /n+i,m+i, P (2)] 

— tbgoyfii ^n)[/n+l,m,p(l) + /n+l,m,p( — 1)] — ^b<?o(<5 + 1, W n ) [/n+l,ro,p+l (0) + /n+l,m,p+l ( — 2)] . (26) 



Here we used the shorthand notations f n (S) = fn(k, S,ui) 
etc. and uj n = uj — ntt in order to simplify the notation. 
Also, 



g (5,u;) = ^J2 eik5aG o( k ^) 



(27) 



are the free propagators in real space, which can be 
calculated analytically] 19 i 36 For any given number n of 
bosons, the needed Green's functions are f n ,m,p(k,S,oj) 
with \S\ < 2, f n ,m(k,S,u}) with \5\ < 3 and f n (k,S,uj) 
with —3 < 6 < 4. Once we know these, we can calculate 
all other / and F generalized Green's functions. 

To solve this set of recurrence equations, we order all 
Green's functions with a given n in a vector V n of di- 
mension In + 1 + 5n(n — l)/2. Then, for any n > 1, 
Eqs. (|14l) - (|26p map into the matrix recurrence equations 



(28) 



V n = a n (k,u)V n -i + /3 n (k,uj)V n+1 



where a n {k,ui) and /3 n (/c,w) are sparse matrices which 
can be read directly from Eqs. (fT4"]) -([2l) )) . The solution, 
for any n > 1, has the general for m 19 i 20 

V n = A n {k,u)V n -x , (29) 

where A n {k,uS) are given by the continued fractions 

a n (k,uj) 



A n {k,uj) 



1 - f3 n (k,uj)A n+1 (k,uj) 



(30) 



starting from a large value N with Apf(k,u) = 0. The 
physical motivation for this choice has been discussed at 
length elsewhere. 19 - Briefly, it is because clouds with too 
many bosons N — ¥ oo are too expensive energetically, and 



I 

therefore very unlikely to be observed. Hence, the propa- 
gators into such states must become vanishingly small for 
a large enough N. In practice one increases N until the 
matrices A n (k,Lu) are converged to within any accuracy 
one chooses. 

Consider now Eq. (|2l?|) for n — 1 . The entries in Vo are 
the various f (k,6,u) = e ikSa G{k,uj) (see Eq. (JTTJl; the 
functions f n ,m an d fn,m,p are defined only for n > 2). As 
a result, once we know the matrix Ai(k, ui), we find 



/i(*,±l,w) = a ± (k,u)G(k,u), 



(31) 



where a ± (fc,a>) are combinations of the appropriate ma- 
trix elements of Ai(k,uj) and e lkSa factors. Using this in 
Eq. (TrjJ gives us the standard solution 



G(k,u) 



1 



to + ir) — efe — ui) 



(32) 



with a self-energy, at this level of MA approximation, of 



S(fc, ui) — — tb [a + (k, ijS) + a (k,u>)^ 



(33) 



In some limiting cases the model has certain symme- 
tries that will insure that many of the generalized Green's 
functions vanish identically. A simple example is the 

tf 

£k = — 2tf cos(fca) — > 0, and so: 



case tf = 0M- Then the free fermion dispersion vanishes, 



go(6, w) -> Ss,o- 



1 



IT) 



(34) 



This simplifies Eqs. (jl4)) - (l26l) significantly and one can 
show that many of the / functions are identically zero. 
For example, only states f 2n , n and f 2n +i,n, /zn+i.n+i sur- 
vive. This becomes obvious if one considers what kinds 
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of boson clouds can be formed if one acts with V repeat- 
edly on any one-fermion state. By continuity, one expects 
that for small tj, these functions will still be very small 
although not precisely zero. This fact could be used to 
speed up the calculation in this limit, by removing these 
kets from the variational set and thus decreasing the size 
of all these matrices. 

In fact, it is known that go(8,ui) decreases expo- 
nentially with increasing distance \8\ for values lu < 
—2dtf, i.e. below the free fermion continuum (c? is the 
dimension) ^ If we are interested in low-energy proper- 
ties at lu < —2dtf, and especially since ui n = lu — ntt 
appears in Eqs. (fT4|) - (|26j) . we see that most of the terms 
on the right are multiplied by exponentially small func- 
tions if \8\ > in their corresponding g (S,uj n ) factor. 
This explains the earlier statement that the approxima- 
tion ignores only exponentially small contributions. One 
can check that going to more extended clouds will bring 
in more terms in Eqs. (|T4l) - ([26|) . but their g (S,Lu n ) fac- 
tors will be even smaller. This argument only becomes 
problematic when the ground-state energy is not too far 
below the free- fermion continuum and Q,/tf —¥ 0. In this 
case, go(S,Lo n ) still decays exponentially with increasing 
8 but very slowly, and much more extended clouds may 
form with significant probability. For the Holstein model, 
this leads to quantitative discrepancies in the MA predic- 
tion for intermediate couplings, where the polaron cloud 
can extend over many sites . 37 ' 38 At strong couplings, the 
Lang-Firsov approach 3 ^ gives a small, one-site cloud po- 
laron, and MA with a one-site cloud restriction becomes 
asymptotically exact. No exact asymptotic solution is 
known for the Edwards model, therefore there is no guar- 
antee that in the limit Q/t/ — > 0, a 3-site cloud will pro- 
vide a good description. However, as we show below, 
MA allows us to decide, with a high level of confidence, 
whether it provides a reasonable description. 

Finally, let us comment on why we chose to allow a 
three-site cloud, as oppose to a one-site or two-site one. 
Consider what would happen if we restricted the varia- 
tional space to single-site clouds only. Then, we can set 
all functions f n , m and f n ,m,p identically to zero. From 
Eq. flT4"l) . we see that the recurrence relation that now 
links together only the f n (8) with 8 — —1,0, 1 does not 
contain any dependence of k in any of its terms, result- 
ing in a self-energy that is independent of k. For exam- 
ple, for tf — this would mean that the polaron is also 
dispersionless. We know that this is not the case due to 
three-site boson processes,— therefore we need to include 
at least such cloud structures in the calculation to cap- 
ture this. If Q is not too small, one can argue that much 
more extended clouds are unlikely for energetical reasons: 
bosons far from the particle cost energy f2 to create yet 
are unlikely to interact with the particle because of the 
separation. Of course, one can increase the size of the 
cloud systematically until convergence is achieved. Such 
an exercise allows one to uncover the physics essential for 
explaining the properties of the polaron, from the impor- 
tance of various terms play. 
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FIG. 3. (Color online) Spectral weights A(k,ui) for k — 
(top), k = ? (middle) and k = n (bottom) from MA (black 
thin lines) and ED (red thick lines). The MA results are 
shifted upwards to ease the comparison. Parameters are 
Q/tb = 2, tf/tb = 0.1, r)/t b = 0.02. 
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FIG. 4. (Color online) Polaron energy E(k) in units of tb 
vs. k, for the parameters indicated. Lines show ED results, 
symbols show MA results (dashed lines are guides to the eye). 



B. MA results for the ID Edwards model 

We can analyze the accuracy of the MA approximation 
for the Edwards model in ID, where accurate numeri- 
cal results have been obtained using variational Exact 
Diagonalization (ED). 2 - We begin with the correlated 
transport regime tt/tb <C 1 (cf. Fig. 1 in Ref. l28r i. 
where polaron motion becomes possible only through 
emission/absorption of bosons. This will be the regime 
of main interest to us when discussing the 2D problem. 

In Fig. [3] we compare spectral weights A{k,Lu) at k = 
0, ? and 7r for f2/t& = 2, tf/tb = 0.1. The agreement 
between MA (thin line, shifted upwards) and ED (thick 
line) results is very satisfactory, especially for the features 
with larger weights. The agreement is of similar quality 
throughout the whole Brillouin zone (not shown). The 
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most prominent features is the low-energy polaron band, 
which disperses very little (its bandwidth is comparable 
to the smallest energy scale in the problem, f/), and has 
a qp weight that varies little with k. 

For a more detailed comparison, we focus on the po- 
laron band, and plot its energy E(k) vs k in Fig. H) As 
expected, for tf — the dispersion corresponds to pure 
2nd NN hopping — 2t2 cos(2fca), where ti is dynamically 
generated through the 3-site, 3-boson processes. For 
tf =/= 0, an additional term — 2t j cos(fca) with a renor- 
malized transfer amplitude tl is also present. Figs. @Ji,b 
show quite good agreement between MA and ED. As ex- 
pected, the agreement is better for the larger 0,/tb = 2 
value, where the probability of more extended clouds is 
reduced. However, even for the smaller fl/tb = 1, MA 
captures the dispersion quite accurately: most of the dif- 
ference to the ED results is an overall shift independent 
of k. This suggests that more extended clouds will not 
further renormalize the effective hopping integrals, only 
lower the overall polaron formation energy. This is rea- 
sonable, because longer range effective hopping terms are 
more complicated to generate and involve many more 
sites than the 3-site, 3-boson process responsible for t^. 
It is therefore probably safe to conclude that MA accu- 
rately describes the polaron's dynamics in the tf/tb <C 1 
region, so long as fi/i/ is not too small so that the spatial 
Green's functions decay exponentially reasonably fast. 

As il/ti, is decreased we expect to move in the regime 
of "strong fluctuations"^ The low-energy part of the 
spectral weight in this regime is shown in Fig. [5] for 
0,/tb = 0.5, tf/tb — 0.04. Because the cost of exciting 
bosons is reduced, we expect to see many more lower- 
energy features than in Fig. [3J and this is indeed the 
case. MA shows reasonable agreement with ED at the 
lower energies, however, some of the higher-energy peaks 
are either missing, or displaced, or combined in a sin- 
gle feature. In order to properly describe these peaks, 





(Color online) Spectral weights A(k,oj) for k — 
= f (middle) and k = it (bottom) from MA and ED. 



FIG. 5. 

(top), k 

The MA results (black thin lines) are shifted upwards 
parameters are fifth = 0.5, tf/tb = 0.04, r)/tb = 0.002 



The 



FIG. 6. (Color online) Spectral weights A(k, uj) for k — 
(top), k = (middle) and k = n/6 (bottom) from MA and 
ED. Again the MA results are shifted upwards. The param- 
eters are Q/t b = 0.5, tf/U = 4, r\/tb = 0.006. 



it is likely that one would need to go to the MA^ or 
higher level approximations, where bosons are allowed to 
exist far from the main polaron cloud (cf. Fig. HJb)). 
Such states are essential in describing the polaron+one 
boson continuum starting at E gs + f V 2 ' 20 and indeed, it 
is at these energies that the disagreements between MA 
and ED become more visible. If, however, the focus is 
on understanding the polaron band and if the polaron 
bandwidth is less than (as is the case we discuss in 
2D, below), then inclusion of these states is not abso- 
lutely necessary: doing so will improve the quantitative 
agreement, of course, but will not change qualitatively 
the polaron dispersion. 

Next we explore the "incoherent" or "diffuse" region 
of the parameter space, where tb/fl ^> 1 while h/tf <C 
1.— Since this implies fi <C tf, this is where the MA 
approximation is expected to be least accurate. In Fig. [6] 
we show comparisons of the spectral function A{k,uj) vs 
uj for small values of k. At k = (upper panel), ED 
shows two sharp peaks, associated with the polaron and 
the second bound state ; 12 ! 20 followed by the polaron+one 
boson continuum at E gs + i7, and then more features 
at higher energies. MA finds the two peaks associated 
with the bound polaron states, shifted to slightly higher 
energies, but the continuum is absent since, as discussed, 
it is not included in the variational space at this MA level. 
As k increases, the spectral weight (equal to the area 
under the peak) in the polaron band decreases extremely 
fast, similar to what happens for Holstein polarons at 
weak coupling ) 11 ! 20 ! 40 ! 41 and most of the weight shifts to 
roughly to = . The lower panels of Fig. [6] illustrate 
this fast decrease in the polaron band spectral weight as 
k increases^ 1 - We also see that even though the polaron 
band predicted by MA is shifted upwards, this shift is 
again not strongly k dependent. The qp weight is also in 
good agreement with the ED results, so MA is still doing 
a reasonable job in predicting the qp weight and effective 
mass. However, at larger fc, the absence of the continuum 
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k/jt 



0.5 
k/ji 



FIG. 7. (Color online) Polaron energy E(k) in units of tb vs. 
fe, for the parameters indicated. Dashed lines show MA re- 
sults, full lines correspond to ED results. The arrows roughly 
mark a discontinuity in the slope of the dispersion. 



within this level of MA approximation is expected to lead 
to an overestimate of the polaron bandwidth*^ 

This is indeed the shown in Fig. where we 

compare the polaron dispersions E(k) vs. fc for 17 = 0.5tf> 
and t f /t b = 2,4,20, so that Q/2t f = 0.125,0.0625 re- 
spectively 0.0125. The agreement at large k worsens as 
il/tf — > 0. This is because in all these cases the band- 
width of the MA E{k) exceeds f2, which is unphysical: 
the polaron band must always stay below the continuum. 
It follows that in this regime, the polaron + one bo- 
son continuum plays the key role in defining the polaron 
bandwidth. To describe it within MA, one needs to allow 
at least one boson to be far away from the polaron, i.e. to 
go to the MAW level. The simpler MA^ ) level approx- 
imation used here is certainly untrustworthy at large k 
in such cases. However, it still provides reasonably good 
description near k = 0, as shown in Figs. [6] and [71 For 
example, the ED polaron dispersion shows a "kink" at 
small k, whose origin we do currently understand. This 
feature is more pronounced as Sl/tf decreases, and is in- 
dicated by the arrows in Fig. [JJ MA also replicates it, 
albeit at somewhat different location. Below this value, 
the two dispersions seem to match quite well (up to an 
overall shift). We conclude that even in this least favor- 
able regime and for this lowest level of MA, we can still 
use it to get reasonable estimates for the polaron effec- 
tive mass and ground-state qp weight. This is verified 
by the data presented in Table U where we compare the 
effective masses (in dimensionless units), 



1 



1 d 2 E(k) 



a 2 tb dk 2 



(35) 



fe=0 



calculated with the two methods. 

To summarize, these comparisons confirm that if the 
MA polaron bandwidth is less than f2 so that the contin- 
uum can be ignored, then this level of MA approximation 
suffices to describe E{k) with good accuracy throughout 
the Brillouin zone. This is generically expected to hold 
so long as tl/tf is not too small, in other words for most 
of the parameter space. In the limit of small Cl/tf, the 
lowest MA approximation suffices for a reasonable de- 





*/ 


m* ED 


m* MA 


1.0 


0.00 


17.71 


18.12 


1.0 


0.01 


15.75 


15.51 


2.0 


0.00 


129.28 


129.98 


2.0 


0.01 


44.68 


44.77 


0.5 


2.00 


1.54 


1.39 


0.5 


4.00 


0.61 


0.47 


0.5 


20.00 


0.052 


0.041 



TABLE I. Comparison of effective masses as predicted by 
ED^ and MA, for several values of and tj where tb = 1. 



scription of GS properties, but one needs to go to higher 
MA levels and enlarge the variational space if one wants 
to have a good description of E(k) in all the Brillouin 
zone. The same holds true if one wants an accurate de- 
scription at energies above the polaron band. Because 
the MA approximations generically satisfy with good ac- 



curacy spectral weight sum rules 



19,20 



one expects even 



the lowest MA level to identify quite correctly where sig- 
nificant spectral weight appears in the spectrum. This is 
indeed the shown in the comparisons provided 

here (regarding Fig. [7j one must remember that there is 
essentially no spectral weight in the qp band once it gets 
close to the continuum). However, in order to capture 
finer details, one needs to work harder by suitably enlarg- 
ing the variational space. In practice, this requires one to 
figure out the corresponding equivalent of Eqs. (|14p - (|2"6")l 
and find an efficient way to solve them. 

In the following, we focus on the polaron dispersion of 
the 2D Edwards model in the limit tf /tb — > for a finite 
0,/tb ratio. In this case, the MA level introduced here — 
suitably generalized to 2D — is sufficient for an accurate 
description of E(k) in the entire Brillouin zone, therefore 
we do not need to go to a higher level. 



C. MA for the 2D square lattice 

To generate the MA equations — within the 3-site cloud 
variational space — for a 2D square lattice, we follow the 
steps outlined in the previous section. The only compli- 
cation is that now the 2-site clouds can be aligned either 
along the x or y directions, so we need two types of 2-site 
cloud generalized Green's functions, 

Fit (k-q,*) = E eJ(k " q)Rl (°l c ^( w )<6j m 6 t— 1 0) 

(36) 

with the associated 



E 



N 



G^cl^l ' », (37) 
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where 77 = (1,0) = x or 77 = (0, 1) = 
is two-dimensional, e.g. i — (i x ,i y ) 
shorthand notation i + 77 = {i x + 1, ij 



y, each site index 
and we used the 
,) if 77 = x etc. As 

I 



in ID, we ask that 1 < m < n — lto avoid overlap with 
the 1-site cloud functions. 

For the 3-site clouds, we define: 



n,m,p V ' Hi 



E e * (k - 



q)R, 



(0|c k G(a;)cUt; 



( |7l— m— p 



(38) 



(k, 5, W ) = - E (k, q,u,) = X; 



ikRi 



^10) 



(39) 



To describe collinear clouds we take 77 = 77', and we 
again need only keep 77 = x, y for the two possible orien- 
tations. For non-collinear clouds we have 4 different dis- 
tinct possibilities which we choose as 77 = ±x = (±1,0) 
and 77' = ±y = (0, ±1), as sketched in Fig. [SJ They are 
inequivalent except when m + p = n, so that there are 
bosons only on the opposite diagonal sites. Again, only 
m > 1, p > 1 are allowed. 

The equations for the various F and / functions are 
generated just as in the ID case, using the restriction that 

I 



I 

boson addition contributions cannot extend the boson 
cloud to more than 3 neighboring sites. This again leads 
to a set of equations for a finite number of n-boson func- 
tions / which depend only on n—1 and n+1 boson / func- 
tions. While straightforward to obtain, these equations 
are very lengthy (e.g. the equation for / n (k, 5, w) has 64 
contributions on the right-hand side), and we do not list 
all of them here. As an example, the relevant equations 
for the 3-site cloud functions are listed below, using again 

the shorthand notations fnJmJ(S) = fn,m,i >(k, S, u) and 

ujn = u) — O: 



=-"»*6 E 9o(6 + S'- V , Un )f<*Zl_ liP (r,)-(n-m-p)t b E 9o(S + S', w„)/^, p (0) 



-VH E 9°(* + S ' + V, w »)/n-?m*-l(-V) ~ t b9o (S + 7,, Un ) E fn+tUl,p(V + S ') 
S' — ±x,±y 8' — ±x,±y 



t b g (S,u n ) E fn+i,L,p( s ') -tbg (8 -r}',u„) E fn+i,L, P +i( s ' ~ v') 



(ViV 1 ) 



5'=±x,±y 



<5'=±x,±y 



(40) 



An analysis of the terms occurring on the right-hand 
sides of these equations allows one to identify the mini- 
mum set of values 5 needed for the various / functions. 
These are shown explicitly for a collinear and one non- 
collinear cloud in Fig. There, the thick lines mark the 
position of the cloud, which is centered at i, while the 
dots show the locations of the fermion. The distances 
from the fermion to the site i define the needed set of 
5 values for these clouds. The other 3-site clouds have 
sets related by appropriate symmetries. The sets for the 
2-site and 1-site clouds are found similarly. 

Once this is done, the solution follows that used in ID. 
All functions with a given n are collected in a vector V n . 
The equations of motions again reduce to matrix recur- 
rence equations V n — a n (k, iv)V n -i +/3 n (k, uj)V n +i which 
are solved in precisely the same way. Of course, the di- 
mension of V n is now significantly increased, dim(V„) = 
31n 2 — 15n — 3, and therefore the various matrices 



r 



A ni a n , f3 n needed are larger than in ID, however they 
are still very manageable. Most results shown below con- 
verged with relative errors less than 10~ 4 if we started 
from An = with N = 9 or less, so that the largest 
vectors' dimension is below 2000. Moreover, their di- 
mensions decrease fast as n decreases, so the solution is 
still very efficient. In the plots shown below, a data point 
typically takes around a minute or less to generate. 

We begin with a thorough analysis of the most in- 
teresting case, when tf = 0. In this case, only the 3- 
site, 3-boson terms already discussed will lead to dy- 
namic generation of a polaron dispersion. We already 
know from the ID case that we expect the generation 
of terms of the type ~ — 2ts[cos(2k x a) + cos(2k y a)] from 
the collinear clouds (on the 2D lattice, this corresponds 
to effective 3rd NN hopping, hence £3). However, be- 
cause of the closed path (Trugman loops) processes that 
are now also possible, we also expect dynamic genera- 



11 



1+Tj 



i-H i 

T1=(1,0),T1'=(0,1) 



l+Tj 



T1=(-1,0),T1'=(0,1) 

i-n 



1+T| 

T1=(1,0),T1=(0,-1) 



i+rj 

n=(-i,o),n'=co,- 



FIG. 8. (Color online) Sketch of our choice of rj, r/' indices 
for the n-boson non-collinear 3-site boson clouds on the 2D 
square lattice. In all cases, m > 1 bosons are at site i — rj, 
p > 1 bosons are at site i + rj' and the remaining n — m — p 
bosons are at site i. 



tion of second NN hopping, leading to terms of the type 



-2t 2 [cos((k x + k y 



cos((k x — k y )a)]. Altogether, 



then, in this case the polaron dispersion should be well 
described by 



£7(k) = E P - 2t 2 [cos((k x + k v 



cos{{k x 



-2ts [cos(2k x a) + cos(2k y a) — 2] 



)a)-2] 
(41) 



where Ep = E(k = 0) is the polaron ground state en- 
ergy. From the ID analysis, we expect that t 2 and £3 
should be quite accurately predicted by MA, while Ep is 
only accurate at fairly high £1 and becomes systematically 
underestimated as f2 decreases. 

In Fig. [10] we plot the polaron dispersion E(k) and qp 
weight Z(k) along lines of high symmetry in the Brillouin 
zone. Interestingly, the two curves have similar profiles, 
however the qp weight changes very little in real terms. 
This is somewhat reminiscent of Holstcin polarons in the 
strong-coupling limit, which also have an almost constant 
qp weight throughout the Brillouin zone. However, in 
that limit their qp weight is exponentially small and the 
effective mass is exponentially large, whereas here the 
qp weight is still considerable, as is the polaron band- 
width. This shows that very different physics gives rise 
to this behavior. Indeed, the small Holstein polaron has 
a small-size cloud with a large average number of bosons. 




FIG. 9. (Color online) Pictorial description of the minimal 
sets of values d needed for a collinear (left) and non-collinear 
(right) 3-site cloud generalized function f n , m ,p{S). Red solid 
lines indicate the boson cloud centered at i, while the blue 
dots mark the possible locations of the fermion. The distances 
from the fermion to the site i define the needed set of S values. 
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FIG. 10. (Color online) Polaron dispersion -E7(k) (upper 
panel) and qp weight Z(k) (lower panel) along the high- 
symmetry directions in the Brillouin zone of the 2D square 
lattice, for tf = 0, tb = = 1. Lines show MA results, the 
symbols are fits to Eq. (|41[) . 



This explains its very large effective mass (due to van- 
ishing overlap of the spatially small polaron cloud), the 
very small qp weight (free fermion contribution to the 
wavefunction is very small), and its weak sensitivity to 
k (states that are nearly localized in real space are "ex- 
tended" in k-space). In contrast, for the Edwards model 
polaron at tf = 0, all the dispersion is due to the exis- 
tence of bosons through the boson-assisted hopping. As 
illustrated in Fig.[TJ the free fermion state mixes with the 
various fermion+bosons states to give rise to the effec- 
tive 2nd and 3rd NN hopping, so the fairly significant qp 
weight throughout the Brillouin zone is not surprising. 
It is worth emphasizing again that the doubling of the 
Brillouin zone is due to the boson-modulated hopping. 
In fact, the resulting dispersion is somewhat reminiscent 
of that of a doped hole in a cuprate layer, although there 
the minimum is at (■£-, ^-), which here is a saddle point. 

The symbols in Fig.QJjlare fits to Eq. (|4Tj) . Specifically, 
we used the MA values for E(0, 0), £7(§ , f ) and £7(0, tt) 
in order to extract Ep,t 2 and £3 from Eq. (|4"T1) . and use 
these to generate the dispersion-fit in the entire Brillouin 
zone. The agreement is very reasonable, backing up our 
assumptions about polaronic physics in this regime. The 
only problem is that £3 <C t 2 ,Ep, and as a result the 
level of confidence in extracting this parameter is not very 
high. For example, if we use £7(0, instead of £7(f , f ) 
as the 3rd point, the value of £3 changes from 0.0023 
to 0.0018 (Ep,t 2 remain unchanged) but the agreement 
between the overall fit and £7(k) is visibly poorer. 

In Fig. [TTJ we study the dependence of t 2 ,t 3 and £7p 
on £l/tb when tf = 0. The full lines show converged MA 
results (except for very small O, see below), whereas the 
symbols show the MA results with the restriction that we 
only allow clouds with up to 3 bosons (A4 = 0). As ex- 
pected, at large f2 the agreement is very good: we know 
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FIG. 11. (Color online) Effective t2/h in (a), tz/tb in (b) and 
Ep/tb in (c), vs. Q/tb for tf — 0. The thick lines give fully 
converged MA results, the dots corresponds to the 3-boson 
solution. The insets show the 3-boson results over a larger 
Q/t range, and the dashed lines show 1/x 5 dependence. See 
text for more details. 
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FIG. 12. Effective hopping integrals £2, £3 and polaron 
ground-state energy Ep vs. the maximum number of allowed 
bosons N. The upper (lower) panels correspond to Q/tb = 0.4 
(Q/tb = 1). In both cases, tf = 0. 



that we need clouds with at least 3 bosons to generate 
the effective hoppings, and because Q is large, it is very 
unlikely to have larger clouds. This is further confirmed 
by the insets, which show that in the limit tb/Q — > 0, 
both effective hoppings scale like t^/Q 5 , as expected for 
the 3-boson, 3-site processes from perturbation theory. 
The fits also confirm that in this limit, t%/t% = 4. This is 
because there is constructive interference in going about 
the Trugman loops clockwise and anticlockwise to gener- 
ate £2, while £3 can only be generated in a unique way. 
As Q decreases below roughly 2£(,, we see that £2, £3 be- 




FIG. 13. Polaron dispersion E(k) (upper panel) along high- 
symmetry cuts in the 2D square-lattice Brillouin zone, for 
tf = 0.1, tb = Q = 1. Lines show MA results, the symbols 
are fits to Eq. (|42|l . The lower panels show the relative error 
e(x) — 1 — x(N)/x(oo) in the effective hopping integrals, as 
well as the polaron ground-state energy, vs. the maximum 
number N of allowed bosons. 



come, at first, larger than the corresponding 3-boson val- 
ues. Indeed, here we have to use larger clouds to achieve 
full convergence, and processes with more than 3 bosons 
will further increase the effective hoppings. Surprisingly, 
for Q/tb < 0.7 or so, t 2 and £3 start to decrease fast. Here 
many-boson processes lead to a decrease of the effective 
hoppings from what the simple 3-boson scenario would 
predict. Convergence of the various quantities in depen- 
dence on the maximum number N of bosons allowed in 
the cloud is shown in Fig. [T2]for Q/tb — 0.4 (upper pan- 
els) and for Q/tb = 1 (lower panels). The data indeed 
confirm that 4 and more boson terms have different ef- 
fects on the effective hoppings for Q < tb and Q > £&. 

It is also clear that in the limit Q/tb — > our results 
are untrustworthy, because we ignore longer loops that 
also contribute to the effective hoppings in this regime. 
For example, just as the 6-step Trugman loop on a 2 x 2 
plaquette contributes to 2nd NN hopping, the Trugman 
loops on 2 x 3 plaquettes will contribute to both 2nd 
and 3rd NN hopping (depending on how the fermion goes 
around the loop) and these contributions will supplement 
the values obtained from the processes included here. Its 
contributions scale asymptotically like t^/Q 10 because 
they involve a 5-boson cloud and 11 hoppings to first 
create and then annihilate all of them. This is to be 
contrasted to t® /Q 5 scaling for the contributions included 
in the 3-site cloud MA. Clearly, once Q ~ tb we cannot 
ignore the contribution of these longer loops. This is 
why it is also pointless to go to larger N to find the fully 
converged values for the Q/tb = 0.4 in Fig. [T3] However, 
a good estimate of the crossover is difficult to obtain from 
such perturbation theory arguments, since the insets in 
Fig. [11] reveal that the asymptotic expressions are only 
valid at much larger Q/% values. 
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FIG. 14. (Color online) The effective hopping integrals and 
the polaron energy, in units of tb, vs. Q,/tb, for tf = O.ltt 
(black solid line) and tf = 0.2tb (red dashed line). The results 
are only shown for SI values where 5 or more boson processes 
become irrelevant. 



A better criterion is to take the value ft j% above which 
convergence is achieved for a cutoff equal to or less than 
6, signaling that 5 or more boson processes are not con- 
tributing much to the polaron wavefunction, and there- 
fore longer loops can be ignored safely. This definition is 
not going to produce a very sharp value since these con- 
tributions change gradually. From Fig. [T2] we see that 
for fl/ti, — 1, the change in going from N = 4 to N = 5 
modifies various quantities by up to about 4%, therefore 
this is likely already in the regime where longer loops are 
not playing an important role. Once 5-boson processes 
become important, we need to include in the variational 
calculation at least the 3x2 loops which will modify both 
ti and t 3 . This is why the numbers shown in Fig. 1111 arc 
likely not valid for small Q. This is an example of how 
a MA approximation can signal its potential problems, 
but also how to fix them (here, extension to at least 5-site 
polaron clouds is needed at lower f2). 

In any case, the boson-modulated hopping is responsi- 
ble for dynamically generating a finite effective (dimen- 
sionless) mass m* = tb/(fa + 2t 3 ) > 14 or so, even though 
the free particle has an infinite mass. 

We next discuss the case of finite tf, in the limit tf < 57 
and fl/tb > 1 where the results of this approximation are 
expected to be valid. The main change in the polaron 
dispersion is that it will also acquire NN contributions, 
so now 

E(k) = E P - 2t f [cos(k x a) + cos(k y a) - 2] 

— 2t 2 [cos((fc ;r + k y )a) + cos((k x — k y )a) — 2] 
-2i 3 [cos(2fc x a) + cos(2k y a) - 2] , (42) 

where t*j is renormalized due to polaron cloud overlap. 

The polaron dispersion is shown in Fig. Q2] for tf — 
O.ltb, f2 = tb- A comparison with Fig. [TU] reveals that 



the small tf has a significant effect on -E(k), especially 
near the edges of the Brillouin zone. The symbols show 
fits to Eq. (|4*2")) . Here, the values for Ep,t* f ,t 2 and t 3 
were extracted from the energies at the special points 
(0,0), (§, |), (0,7r), and (n,ir). Using these parameters 
in Eq. (|4"2l leads to good agreement with the MA re- 
sults (thick lines). We find that t* f /t b = 0.040, t 2 /t b = 
0.0090, t 3 /t b = 0.0026. Interestingly, while t* f < tf, as 
expected in polaronic physics, we see that in the presence 
of a finite tf, both t 2 and t 3 are larger than when tf = 0. 
This is because the effective hopping integrals generated 
by the 3-boson processes discussed so far are here supple- 
mented by the usual longer-range polaron hopping known 
to occur in the intermediate-to-strong electron-phonon 
coupling limiti 40 ' 41 Just as for tf ~ 0, the qp weight 
changes little throughout the Brillouin zone; it is close 
to 0.3 everywhere. The lower panels in Fig. [T3]show the 
convergence of the various effective hopping amplitudes 
and of the ground-state energy as the maximum number 
of bosons N increases. Specifically, we plotted the rela- 
tive errors which reveal that 5 or more boson processes 
contribute less than 4% to the various quantities. 

The dependence of t*j,t 2 , t 3 and Ep on il/t is shown 
in Fig. [14] for a two values of tf. The data is only dis- 
played over the range where 5-boson processes contribute 
less than 1% to the various parameters, so that longer 
Trugman loops can be safely ignored. As expected, t*j 
increases towards tf as increases, because this leads to 
fewer bosons in the polaronic cloud, i.e. less "dressing" of 
the quasiparticle. On the other hand, t 2 and £3 decrease 
with increasing Q, as this makes the intermediary many- 
boson states less likely. As already discussed, their values 
increase with increasing tf, for a fixed value of Q. Finally, 
we note that Ep is well below the free particle continuum 
starting at —4tf. This, together with the fact that longer 
loops are irrelevant, guarantees that the approximation 
must be quantitatively accurate in this regime. 



IV. SUMMARY AND DISCUSSIONS 

The main goal of this work is to demonstrate how the 
MA approximations can be generalized to study polaron 
formation in models with boson-affected fermion hop- 
ping. Unlike for simpler local fermion-boson couplings as 
in the Holstein polaron model, where the strong-coupling 
Lang-Firsov solution is known and can guide the choice 
for the maximum extension of the polaron cloud, here this 
solution is not available. As illustrated for the Edwards 
fermion-boson model, one now needs to use physical in- 
tuition to make a reasonable choice. Of course, one can 
always systematically increase the variational space and 
check that the initial guess was indeed reasonable. 

Because we are interested primarily in the low-energy 
polaron physics, we used the MA' - 1 level of approxima- 
tion which describes only the polaron cloud and does not 
allow for far-flung bosonic excitations. Then, the only 
free "parameter" is the spatial size of the polaronic cloud. 
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Simple arguments regarding the processes illustrated in 
Fig. [1] show that at least 3-site clouds need to be al- 
lowed, and therefore we built the approximation for this 
case. Careful consideration of the terms ignored gives 
us intuition about when the approximation is expected 
to be accurate, and when and in which way it becomes 
problematic. In ID, this was indeed verified successfully 
against available exact numerical results. 

We then extended the calculation to 2D, where no re- 
sults were available for this model until now, and demon- 
strated that the closed Trugman loops play the key role 
in determining the effective mass of the quasiparticle in 
the limit tf = 0. In this regime, the results for the 3-site 
MA calculation are trustworthy for SI > f&; for smaller f2 
values, one needs to increase the allowed size of the po- 
laron cloud since longer Trugman loops are also becoming 
important. We emphasize that the MA approximation is 
not wrong in the low Q limit; what failed is our guess 
about the relevant size of the polaron cloud. If this is 
increased from 3 to more sites, the approximation will 
become accurate again. 

The important role played by Trugman loops raises a 
very important question regarding the motion of a par- 
ticle through an AFM background (which, as discussed, 
the Edwards fermion-boson model partly mimics). For 
t — J models, it has been argued that the interaction of 
the hole with spin- waves is well described within the self- 
consistent Born approximation. This approximation in- 
cludes only non-crossing diagrams, i.e. processes where 
the bosons are absorbed in inverse order to the one in 
which they have been emitted by the particle. This is 
because it is generally expected that the particle needs 
to retrace its path to "heal" the string of defects it cre- 
ated when it reshuffled the spins. For a Neel AFM, the 



self-consistent Born approximation therefore predicts an 
infinitely heavy quasiparticle. In the presence of spin- 
fluctuations, the magnons can disperse and this gives rise 
to a finite quasiparticle mass. 

What is shown here is that the quasiparticle can ac- 
quire a finite spin mass even in the absence of spin fluctu- 
ations, by going (almost) twice around closed loops, first 
creating a string of defects and then healing it. Note 
that in diagrammatic terms, this would correspond to 
maximally crossed diagrams, since here the bosons are 
absorbed in the same order in which they were emit- 
ted. Such processes are obviously not included in the 
self-consistent Born approximation. 

Of course, one might argue that spin fluctuations are 
relevant for holes in cuprates, in other words tf may be 
considerable and therefore dominate E{k), as we found it 
to be the case for the Edwards model for tf > 0. In other 
words, that the Trugman loops' contributions, although 
finite, may be quantitatively insignificant. However, the 
interesting thing is that these closed loop processes give 
rise to the only contributions in E{k) which are consistent 
with the doubling of the unit cell. ARPES on the parent 
insulators, i.e. for a single quasiparticle introduced in a 
Cu02 layer, clearly exhibits this doubling of the Brillouin 
zone^ This suggests that maybe this problem should be 
revisited using MA-type approximations. 
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